#この枠をクリックしてアクティブにしてから,Shiftを押しながらEnterを押すと,枠内のコードが実行されます.以下同じです.
# '#'記号のあとの内容を「コメントアウト」と呼び,実行時に無視されます.
import numpy as np #数値計算ライブラリ
from math import sin, cos #低機能だが計算が速い三角関数
from scipy.integrate import odeint #常微分方程式ライブラリ
import matplotlib.pyplot as plt #グラフ作成ライブラリ
from matplotlib.animation import FuncAnimation #アニメーションライブラリ
from matplotlib import rc #グラフ調整ライブラリ
from IPython.display import HTML #ユーザーインターフェース拡張ライブラリ
M = 2/3 #台車の質量
m = 1/3 #振り子の質量
l = 1 #振り子の長さ
def Solve(xx0, force):
### 運動方程式を1階化した微分方程式の定義
def eom(xx,t):
#注 python は配列の添字が0から
x = xx[0] #台車の変位
dx = xx[1] #台車の速度
th = xx[2] #振り子の角度
dth = xx[3] #振り子の角速度
ft = force(xx) #制御力
A = np.array( #行列
[
[M+m, m*l*cos(th)],
[np.cos(th), l]
]
)
bb = np.array( #右辺のベクトル
[m*l*(dth**2)*sin(th)+ft, 9.8*sin(th)]
)
A_inv = np.linalg.inv(A) #逆行列
hh = np.dot(A_inv, bb) #逆行列とベクトルの積
dxx = np.array( #1階微分のベクトル
[dx, hh[0], dth, hh[1]]
)
return dxx
###ここまで
### 時間軸を表す等差数列
n = 200 #時刻の数
ts = np.linspace(0, 25, n) #0秒から25秒までn等分
### 差分解法(数値積分)
xxs = odeint(eom, xx0, ts)
return (ts,xxs)
def force(xx):
return 0
xx0 = np.array([0, 0, 0.5, 0]) #時刻t=0 のx, dx, th, dth
ts, xxs = Solve(xx0, force)
ベクトル xx := (x, dx, th, dth) を **状態ベクトル** という.
print(xxs) #各行が状態ベクトル,下に向かって時間が進行
def Plot_Wave(ts, xxs):
plt.plot(ts,xxs[:,0],label='x') #台車の変位
plt.plot(ts,xxs[:,1],label='dx/dt') #台車の速度
plt.plot(ts,xxs[:,2],label='th') #振り子の角度
plt.plot(ts,xxs[:,3],label='d(th)/dt') #振り子の角速度
plt.xlabel('t')
plt.ylabel('States')
plt.legend()
Plot_Wave(ts, xxs)
def Animate(xxs, title='xxxxxxX'):
### アニメーション用のグラフ用紙を用意する
fig, ax = plt.subplots(figsize=(8,3)) #グラフ用紙を作る
plt.close() #ひとまず表示OFF
ax.set_xlim(-4,4) #グラフの縦軸の範囲
ax.set_ylim(-1.5,1.5) #グラフの横軸の範囲
ax.grid() #グリッドon
ax.set_xlabel('X') #横軸のラベル
ax.set_ylabel('Y') #縦軸のラベル
ax.set_title(title) #タイトル(学籍番号)
line1, = ax.plot([], [], lw=2) #空の描画
### i 行目の状態ベクトルで描画データを更新する関数
def update_anim(i):
x = xxs[i,0] #さっき計算したxxsの1列目(台車変位,添字は0)
th = xxs[i,2] #さっき計算したxxsの3列目(振子角度,添字は2)
#振り子支点の位置ベクトル
XM = np.array([x, 0])
#振り子先端の位置ベクトル
Xm = XM + l*np.array([sin(th), cos(th)])
#振り子の線分を描画
line1.set_data([XM[0],Xm[0]],[XM[1],Xm[1]]) #座標データの更新
### アニメーションデータの作成
n = len(xxs) #データの行数(時間きざみの総数)
anim = FuncAnimation(fig, update_anim, interval=100, frames=n)
### アニメーションの出力形式の設定
rc('animation', html='jshtml')
return anim
Animate(xxs)
forceとxx0を,def force(xx):
return 20*xx[2]
xx0 = np.array([0, -3, 0.5, 1.2])
のように書き換えて実行せよ.
def force(xx):
return 0
xx0 = np.array([0, 0, 0.5, 0])
ts, xxs = Solve(xx0, force)
Animate(xxs)
forceとxx0を,def force(xx):
return 20*xx[2] + 2*xx[3]
xx0 = np.array([0, -3, 0.3, 1.2])
のように書き換えて実行せよ.
def force(xx):
return 0
xx0 = np.array([0, 0, 0.5, 0])
ts, xxs = Solve(xx0, force)
Animate(xxs)
次のコードセルのKとLの数値を自分で書き換え(必要ならxx0も),次の2種類の運動を実現せよ.
(1) 振動しながら立位に向う運動.
(2) 振動しないで立位に向う運動.
def force(xx):
K = 20
L = 2
return K*xx[2] + L*xx[3]
xx0 = np.array([0, -3, 0.3, 1.2])
ts, xxs = Solve(xx0, force)
Animate(xxs)